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Abstract: Gravity surveys are an important research topic in geophysics and geodynamics. 
This paper investigates a method for high accuracy large scale gravity anomaly data 
reconstruction. Based on the airborne gravimetry technology, a flight test was carried out 
in China with the strap-down airborne gravimeter (SGA-WZ) developed by the Laboratory 
of Inertial Technology of the National University of Defense Technology. Taking into 
account the sparsity of airborne gravimetry by the discrete Fourier transform (DFT), this 
paper proposes a method for gravity anomaly data reconstruction using the theory of 
compressed sensing (CS). The gravity anomaly data reconstruction is an ill-posed inverse 
problem, which can be transformed into a sparse optimization problem. This paper uses the 
zero-norm as the objective function and presents a greedy algorithm called Orthogonal 
Matching Pursuit (OMP) to solve the corresponding minimization problem. The test results 
have revealed that the compressed sampling rate is approximately 14%, the standard 
deviation of the reconstruction error by OMP is 0.03 mGal and the signal-to-noise ratio 
(SNR) is 56.48 dB. In contrast, the standard deviation of the reconstruction error by the 
existing nearest-interpolation method (NIPM) is 0.15 mGal and the SNR is 42.29 dB. 
These results have shown that the OMP algorithm can reconstruct the gravity anomaly data 
with higher accuracy and fewer measurements. 
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1. Introduction 

The Earth's gravity field is a fundamental physical field, which reflects the distribution, motion and 
variety of the Earth's interior matter. The gravity has connections with all the physical events on Earth 
and in its near space, and thus provides the basic information to support research on many subjects. 
Gravity surveys can support fundamental geophysical investigations, which are beneficial to determine 
the density of the Earth's interior matter and help explain many physical phenomena of the Earth. 
Meanwhile, gravity surveys are also significant in the exploitation of mineral resources and modern 
military science, etc. 

Airborne gravimetry is a method of determining the Earth's gravity by using instruments on board 
an aircraft such as accelerometers, global navigation satellite systems (GNSS), altimeters, and attitude 
sensors [1,2]. Performing gravity surveys from an aircraft is superior to point- wise terrestrial 
gravimetry in terms of both economy and efficiency. It also offers the opportunity to measure gravity 
over special terrains which are difficult to access and in the areas with mixed land and ocean [3,4]. 
Airborne gravimetry can be applied to rapidly get extensively well-distributed information about the 
Earth's gravity field, so it is an important new geodetic technology and an essential part of airborne 
geophysical exploration and also is supplemental to terrestrial gravimetry and shipborne gravimetry. 

Airborne gravimetry is essentially a discrete digital sampling method. The theoretical foundation of 
discrete digital sampling on continuous-time band-limited signals was developed by Nyquist and 
Shannon [5,6]. Their results demonstrated that signals can be recovered from a set of uniformly spaced 
samples taken at the so-called Nyquist rate of twice the highest frequency present in the signal of 
interest. However, because of the restrictions posed by national boundaries, economic cost and 
database size, airborne gravimetry is in practice a sub-Nyquist sampling method. Consequently, there 
is a question on whether the gravity data can be recovered with a new framework for signal acquisition 
and sensor design that enables a potentially large reduction in the sampling and computation costs for 
sensing signals. More recently, Candes, Romberg, Tao and Donoho showed that a signal having a sparse 
representation can be reconstructed from a small set of linear, non-adaptive measurements [7-11]. This 
result suggests that it may be possible to sense sparse signals by taking far fewer flight measurements, 
and hence such as a method is named Compressed Sensing (CS). 

Based on airborne gravimetry technology, a flight test was carried out in China with the 
strap-down airborne gravimeter (SGA-WZ) developed by the Laboratory of Inertial Technology of the 
National University of Defense Technology [12-14]. This paper investigates a method for large scale 
gravity anomaly data reconstruction. Taking into account the sparsity of airborne gravimetry by DFT, 
this paper firstly proposes a method for the reconstruction of gravity anomaly data using the CS theory. 
The gravity anomaly data reconstruction is an ill-posed inverse problem, which can be transformed 
into a sparse optimization problem. This paper uses the zero-norm as the objective function and 
presents the OMP algorithm to solve the corresponding minimization problem [15]. These results have 
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shown that the OMP algorithm can reconstruct the gravity anomaly data with higher accuracy and 
fewer measurements than existing methods 

2. Basic Principle and Mathematical Model 

Airborne gravimetry can be classified as airborne scalar gravimetry, airborne vector gravimetry and 
airborne gradient gravimetry. Our work in this paper studies airborne scalar gravimetry 
(we study vertical component of the vector only). 

2.7. Airborne Gravimetry Technology 

An object's gravity is the composition of forces of gravitation caused by the Earth and other 
celestial bodies and the inertial centrifugal force caused by the Earth's rotation. The non-uniform 
distribution of the density of the Earth's interior matter makes the gravity vary with the position. In 
gravity prospecting, the gravity variations caused by the non-uniform density distribution of the 
Earth's interior rocks and minerals are called gravity anomalies. In fact, the airborne gravity 
measurements include two parts: the gravity anomaly (denoted as the free air anomaly here) and the 
normal gravity (which is the reference gravity field of a conventional ellipsoid). Therefore, the gravity 
anomaly can be expressed as: 

^g = g-r (1) 

where Ag is the gravity anomaly, g is the gravity measurement and y is the normal gravity. 

The principle of strap-down airborne scalar gravimetry is based on Newton's equation of motion in 
the gravitational field of the Earth, utilizing the principle of relative gravity measurement. First, we 
must use a terrestrial gravimeter to connect national gravity and a point on the parking apron to obtain 
its gravity and take it as the gravity reference point. Before the aircraft takes off, static gravimeter data 
must be recorded so that gravity observations in the air can be associated with the gravity reference 
point on the parking apron. Then, the strap-down airborne scalar gravimetry model can be written as: 

Ag = - {fo -fD+8b) + - 7 (2) 

Obviously, Equation (2) is an equivalent form of Equation (1), where Ag is the gravity anomaly to 
be determined, is the vertical component of the vehicle acceleration obtained from GPS, is the 
vertical component of the specific force measured by the accelerometers of an inertial measuring unit, 
fl is the vertical component of the specific force obtained from the static data on the parking apron, 
called base reading, g^ is the gravity of reference point on the parking apron and iSa^ includes all 
kinds of the error correction [16]. 

2.2. Mathematical Model of CS 

Compressed sensing (CS) is an exciting, rapidly growing field that has attracted considerable 
attention in fields as diverse as electrical engineering, applied mathematics, statistics, sensor 
technology and computer science. CS offers a framework for simultaneous sensing and compression of 
finite dimensional signals. Quite surprisingly, it predicts that sparse high-dimensional signals can be 
recovered from highly incomplete measurements by using efficient algorithms. CS also holds promise 
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for increasing resolution by exploiting the signal structure. Especially, reducing the sampling rate or 
increasing resolution in airborne gravimetry can improve survey efficiency, increase data transfer rate 
and improve data quality. 

Let z be an unknown original gravity anomaly in (which can represent a 1-D or 2-D gravity 
anomaly of interest). Suppose that we have m linear measurements of z with the form: 

J.=(^'Z) + v., / = l,2,---,m (3) 

where (•,•) denotes the usual inner product, j = [ji,j2'*"'Jm]^ the sampling gravity anomaly, 

0 = [<z!(,^2'* * '^^mf ^ ^ is known and called the compressed sensing matrix and v = [v^^v^,- • sv^]^ g R"^ 
is the noise. Standard reconstruction methods require at least n samples. Suppose we know a priori 
that z is compressible or has a sparse representation in a transform domain, described by the matrix 
W G . The sparse representation of z can be written as: 

x = 'Fz (4) 

where x is the sparse coefficient matrix and ^ is the sparse transform matrix. In this case, if vectors 
^. are well chosen, then the number of measurements m can be dramatically smaller than the size n 
usually considered necessary [11,17]. 

When W is invertible, CS exploits the sampling gravity anomaly data reconstruction by solving a 
problem of the form: 

minimize s.t.y = 0W~^x (5) 

where ||x||q =|{/:x- ^0}| is the zero-norm of x and denotes the number of nonzero elements within x . 
Suppose X has at most k nonzero elements, i.e., \\x\\q <^, we can consider that z is ^-sparse in the 
transform domain W . 

2.3. OMP Algorithm 

Let A = OW~^ G IT"'' ; Equation (5) can be rewritten as: 

minimize ||x||q s.t. y = Ax (6) 

Equation (6) is the zero-norm sparse optimization problem and can be solved with greedy 
algorithms. Greedy algorithms rely on iterative approximation of the signal coefficients and support, 
either by iteratively identifying the support of the signal until a convergence criterion is met, or 
alternatively, by obtaining an improved estimate of the sparse signal at each iteration that attempts to 
account for the mismatch to the measured data. Greedy algorithms are classified as greedy pursuits and 
thresholding type algorithms. Greedy pursuit algorithms include Matching Pursuit [18], OMP [19], 
Conjugate Gradient Pursuit [20], Stagewise Orthogonal Matching Pursuit [21] and Regularized 
OMP [22,23], etc. 

OMP is a greedy pursuit algorithm. Given that z is A:-sparse in the transform domain ^ .i.e., x has 
at most k nonzero elements and A = =\a^,a.^, " whose columns denoted by a.^K^ , the 
vector J = Ar is a linear combination of k columns from the matrix A . In the language of sparse 
approximation, we say that 3; has an ^-term representation over the dictionary A . Therefore, OMP can 

be used for recovering sparse signals. 
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To identify the signal x , we need to determine which columns of A participate in the vector y . 
The idea behind the OMP algorithm is to pick columns in a greedy fashion. At each iteration, we 
choose the column of A that is most strongly correlated with the remaining part of y . Then we 
subtract off its contribution to 3; and iterate on the residual. One hopes that, after k iterations, the 
OMP algorithm will have identified the correct set of columns. The major advantages of OMP are its 
speed and its ease of implementation. The procedure of OMP can be carried out by means of the basic 
steps as follows [15]: 

Given the matrix A , the vector y and the sparsity level k of the signal x : 

Step 1 : Initialize the residual = 3; , the index set 4) = 0 ^i^d the iteration counter t^l. 
Step 2: Find the index that solves the easy optimization problem: 

l,=argmax.^i...J(r,_i,a.)| (7) 
If the maximum occurs for multiple indices, break the tie deterministically. 

Step 3: Augment the index set and the matrix of chosen atoms: =A^_^\j{X^} and ^ ^^t-i ^ • 
We use the convention that /2q is an empty matrix. 

Step 4: Solve a least squares problem to obtain a new signal estimate: 

s^=argrmn^\\y-n4^=(nfnynfy (8) 
Step 5: Calculate the new residual: 

r, = y-^s^ (9) 

It is important to note that the residual is always orthogonal to the columns of . Provided that 
the residual r^_^ is nonzero, the algorithm selects a new atom at iteration t and the matrix i?^ has full 
column rank. For such a case, the solution s^ to the least squares problem in Step 4 is unique. 

Step 6: Increment t , and return to Step 2 if t<m. 

Step 7: The estimate x for the signal x has nonzero indices at the components listed in . The value 

of the estimate x in component equals the rth component of s^ . 
Step 8: The estimate for the unknown signal z is: z = W^x . 

3. Test Results and Discussion 

The strap-down airborne gravimeter and the flight test results are presented in this section. By 
comparing NIPM for the gravity anomaly data reconstruction, the superiority of the OMP algorithm 
will be shown in this section. 

3.1, Flight Test and Data Preprocessing 

The strap-down airborne scalar gravimeter called SGA-WZ mentioned in this paper is the first 
system with this type in China. It was developed by the Laboratory of Inertial Technology of the 
National University of Defense Technology [13]. This system consists of a high-performance 
strap-down inertial navigation system (SINS), a GNSS receiver, an anti- vibration system, a data logger 
and post-processing software. The major advantage of this system is its reliability and robustness in 
operations. A photograph of the system is shown in Figure 1. 
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Figure 1. Appearance of the SGA-WZ. 




The flight tests were carried out in Shandong Province using SGA-WZ from April 2010 to May 
2010. The hardware was installed in the aircraft six days before the flight tests. Two GNSS ground 
stations were located near the airport where the aircraft took off and landed. GNSS receivers installed 
on the ground and aircraft can be used to determine the vehicle position, velocity, and acceleration. 
The strap-down airborne scalar gravimeter onboard a Cessna 208 aircraft was used to collect the data. 
Figure 2 shows the Cessna 208 aircraft, which is a fixed-wing small aircraft. 

Figure 2. Cessna 208 fixed- wing small aircraft. 



Sensors 2014, 14 



5432 



The pilots controlled the aircraft with the autopilot, and the test was implemented in days with good 
weather to minimize the effects of air turbulence. The average flight altitude was approximately 400 m 
above sea level with a fluctuation of 20 m. The average speed during the flight was 60 m/s. The 
sampling rate of raw SINS readings was 100 Hz and 2 Hz for the GNSS sampling rate. After being 
installed in the aircraft, SGA-WZ worked all day for over a month. In the whole flight test campaign 
there were eight flights. The first and second flights were repeated lines that flew along the same 
trajectory to test the repeatability of the system, and the valid length of each repeated profile was about 
100 km. The other six flights left were grid flights consisting of three flights of survey lines and three 
flights of control lines. The spacing between survey lines was about 2 km and the spacing between 
control lines was about 9 km. Figure 3 shows the grid flight lines used in the test [16]. 

Figure 3. Grid flight lines. 



0.8 



0 

■D 
3 



0.6^ 



05 
_l 
I 

CD 
> 

0 

DC 



0.4- 



0.2 



0.2 0.4 0.6 

Relative-Longitude (°) 



0.8 



As previously mentioned, the gravity anomaly was determined from the difference between the 
specific force and the vehicle acceleration using Equation (2). Because of the phugoid modes of the 
aircraft and the atmospheric turbulence, there were many high frequency noises in the 
gravity anomaly. The noises should be eliminated by the FIR low-pass filter. When the velocity of 
aircraft is fixed, the longer the filter length, the more high frequency noises can be eliminated, but the 
spatial resolution will be lowered. The test results have revealed that the internal accord precision is 
1.50 mGal with a spatial resolution of 6 km. Figure 4 shows the gravity anomaly of line 1 in the 
first flight before low-pass filtering. Figure 5 shows the gravity anomaly after 160 s cutoff 
low-pass filtering. 
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Figure 4. Gravity anomaly before low-pass filtering. 
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3.2, Spar sity Analysis 

The first and second flights were repeated lines surveys denoted by F401 and F402, respectively. 
Each of them consisted of six lines. Figure 6 shows the gravity anomaly curves of F401 and Figure 7 
shows the gravity anomaly curves of F402. In the two flights, the repeated lines were surveyed along 
an east-west direction. 
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Figure 6. Gravity anomaly curves of F401. 
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Figure 7. Gravity anomaly curves of F402. 
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Based on the theory of CS, this study analyzes the sparsity of airborne gravimetry by DFT. Figure 8 
shows the amplitude spectrum of all lines in F401 and Figure 9 shows the amplitude spectrum of all 
lines in F402. From the two figures, we can see that the gravity anomaly is mostly distributed in the 
low-frequency domain and it has few high-frequency components. The results thus reveal that the 
gravity anomaly is sparse in the DFT domain, which satisfies the precondition for application of CS. 
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Figure 8. Amplitude spectrum curves of F401. 
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Figure 9. Amplitude spectrum curves of F402. 
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3.3. OMP Reconstruction Results 

To evaluate the OMP algorithm performance and the gravity anomaly data reconstruction precision, 
we now define the performance indices with standard deviation and SNR. SNR can be calculated by: 

SNR = 201og,o(||4/||z-4) (10) 

In this study, the compressed measurement number was set as m = 432, i.e., the compressed 
sampling rate was approximately 14%. Figure 10 shows the reconstruction result of line 1 in F401 and 
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the reconstruction error is shown in Figure 11. In Figure 10, the red dashed line denotes the 
reconstruction gravity anomaly and the blue solid one denotes the original gravity anomaly. The two 
figures show that the OMP results approximate the original gravity anomaly very well and the 
reconstruction error is very small. However, we also note that a boundary effect exists in the result and 
this problem should be solved in the future. 

Figure 10. Reconstruction result of OMP in F401. 
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Figure 11. Reconstruction error of OMP in F401. 
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For comparison purposes, we also made the gravity anomaly data reconstruction with NIPM, which 
is one of the commonly used methodologies for data reconstruction. Figure 12 shows the 
reconstruction result of line 1 in F401 and the reconstruction error is shown in Figure 13. The 
comparison of the algorithm performance between OMP and NIPM is shown in Table 1 . It shows that 
the each performance index from OMP is superior to the corresponding one from NIPM. The standard 
deviation of the reconstruction error and SNR of OMP are 0.03 mGal and 56.48 dB, respectively. 
Correspondingly, the standard deviation of reconstruction error and SNR of NIPM are 0.15 mGal and 
42.29 dB, respectively. 



Figure 12. Reconstruction result of NIPM in F401. 
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Table 1. The comparison of the algorithm performance for F401. 



Number 


Points 


Standard Deviation (mGal) 


SNR (dB) 


OMP 


NIPM 


OMP 


NIPM 


Line 1 


3050 


0.03 


0.15 


56.91 


42.06 


Line 2 


3150 


0.03 


0.16 


55.76 


41.03 


Line 3 


2979 


0.02 


0.14 


58.15 


42.81 


Line 4 


3120 


0.03 


0.14 


56.65 


43.00 


Line 5 


2964 


0.03 


0.16 


57.08 


42.23 


Line 6 


3170 


0.04 


0.15 


54.31 


42.60 


Average 


3072 


0.03 


0.15 


56.48 


42.29 



4. Conclusions 

This paper has presented a new investigation on a method for high accuracy large scale gravity 
anomaly data reconstruction. Based on the airborne gravimetry technology, a flight test was carried out 
in China using a custom designed strap-down airborne gravimeter. Taking into account the sparsity of 
airborne gravimetry by DFT, this paper proposes a method for gravity anomaly data reconstruction 
using CS theory. The test results have revealed that the compressed sampling rate is approximately 
14%, the standard deviation of the reconstruction error by OMP is 0.03 mGal and SNR is 56.48 dB. In 
contrast, the standard deviation of the reconstruction error by NIPM is 0.15 mGal and SNR is 42.29 dB. 
These results have shown that OMP can reconstruct the gravity anomaly data with higher accuracy and 
fewer measurements. In future investigations, the following considerations should be taken into 
account: 

(1) The boundary effect exists in the result of gravity anomaly data reconstruction with OMP. 
Although we can ignore the boundary data, that will result in a waste of data and increase the 
survey costs. 

(2) The discrete Fourier transform was effective for the 1-D gravity anomaly data sparse transform 
in this study. In consideration of the 2-D gravity anomaly data reconstruction, future work can 
pay attention to the Curvelet transform. 
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